researchSprint
2025-08-21
Last updated: 2025-09-23
Checks: 7 0
Knit directory: methyl_nano_cf/
This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20250606) was run prior to running
the code in the R Markdown file. Setting a seed ensures that any results
that rely on randomness, e.g. subsampling or permutations, are
reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 1e62983. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for
the analysis have been committed to Git prior to generating the results
(you can use wflow_publish or
wflow_git_commit). workflowr only checks the R Markdown
file, but you know if there are other scripts or data files that it
depends on. Below is the status of the Git repository when the results
were generated:
Ignored files:
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: analysis/figure/
Ignored: renv.lock
Ignored: renv/
Untracked files:
Untracked: analysis/filteringReads.Rmd
Unstaged changes:
Modified: .Rprofile
Modified: .gitignore
Modified: analysis/alphaAnalysis.Rmd
Modified: analysis/checkAlphaAndBetaMatch.Rmd
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were
made to the R Markdown (analysis/researchSprint.Rmd) and
HTML (docs/researchSprint.html) files. If you’ve configured
a remote Git repository (see ?wflow_git_remote), click on
the hyperlinks in the table below to view the files as they were in that
past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | 9227458 | caitlinpage | 2025-08-28 | moved bins, genes, and cpg features analysis to researchSprintB to make file less massive |
| Rmd | 79d6a02 | caitlinpage | 2025-08-28 | add description, tidy up code |
| Rmd | ded83b0 | Caitlin Page | 2025-08-25 | downloaded from research sprint server |
Description
- 2 day research sprint on long read data from the 21-22 August 2025
- Team leader: Dineika
- Interests: fragmentomics, methylation, deconvolution, slide
- Team members: Dineika, Max, Jovana, Gunjan, Caitlin, Therese, Song, Kevin
- Team name: Goldilocks
What was done
- beta and pat files from wgbs_tools (Max)
- tumour fraction estimates - ichorCNA (Song)
- Fragle (fragment based deconv - short read designed) (Dineika)
- end motif analysis (Gunjan)
- processing data, variance analysis, genome binning (Jovana)
- correlate tumour fraction and proportion of meth/unmeth sites at luad genes? (Therese - I helped)
what did I do?
- some basic plots
- binning the genome
- plots of luad genes
- cpg islands, gene promoters plots
- tried and failed to make beta files using the alpha, then made again from bedgraphs
summary of sample information
- bc10 and bc11 have the highest tumour fraction
- bc8 and bc9 have the lowest tumour fraction - these 2 are grouped together as low - everything else is high
my analysis code
library(biomaRt)
library(EnsDb.Hsapiens.v86)Loading required package: ensembldbLoading required package: BiocGenerics
Attaching package: 'BiocGenerics'The following objects are masked from 'package:stats':
IQR, mad, sd, var, xtabsThe following objects are masked from 'package:base':
anyDuplicated, aperm, append, as.data.frame, basename, cbind,
colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
Position, rank, rbind, Reduce, rownames, sapply, saveRDS, setdiff,
table, tapply, union, unique, unsplit, which.max, which.minLoading required package: GenomicRangesLoading required package: stats4Loading required package: S4Vectors
Attaching package: 'S4Vectors'The following object is masked from 'package:utils':
findMatchesThe following objects are masked from 'package:base':
expand.grid, I, unnameLoading required package: IRangesLoading required package: GenomeInfoDbLoading required package: GenomicFeaturesLoading required package: AnnotationDbiLoading required package: BiobaseWelcome to Bioconductor
Vignettes contain introductory material; view with
'browseVignettes()'. To cite Bioconductor, see
'citation("Biobase")', and for packages 'citation("pkgname")'.Warning: replacing previous import 'S4Arrays::makeNindexFromArrayViewport' by
'DelayedArray::makeNindexFromArrayViewport' when loading 'SummarizedExperiment'Loading required package: AnnotationFilter
Attaching package: 'ensembldb'The following object is masked from 'package:stats':
filterlibrary(plyranges)
Attaching package: 'plyranges'The following objects are masked from 'package:ensembldb':
filter, selectThe following object is masked from 'package:AnnotationDbi':
selectThe following object is masked from 'package:IRanges':
sliceThe following object is masked from 'package:biomaRt':
selectThe following object is masked from 'package:stats':
filterlibrary(tidyr)
Attaching package: 'tidyr'The following object is masked from 'package:S4Vectors':
expandlibrary(ggplot2)
library(BSgenome.Hsapiens.UCSC.hg38)Loading required package: BSgenomeLoading required package: BiostringsLoading required package: XVector
Attaching package: 'Biostrings'The following object is masked from 'package:base':
strsplitLoading required package: BiocIOLoading required package: rtracklayer
Attaching package: 'rtracklayer'The following object is masked from 'package:BiocIO':
FileForFormatlibrary(stringr)
library(data.table)
Attaching package: 'data.table'The following object is masked from 'package:plyranges':
betweenThe following object is masked from 'package:GenomicRanges':
shiftThe following object is masked from 'package:IRanges':
shiftThe following objects are masked from 'package:S4Vectors':
first, secondlibrary(dplyr)
Attaching package: 'dplyr'The following objects are masked from 'package:data.table':
between, first, lastThe following objects are masked from 'package:Biostrings':
collapse, intersect, setdiff, setequal, unionThe following object is masked from 'package:XVector':
sliceThe following objects are masked from 'package:plyranges':
between, n, n_distinctThe following objects are masked from 'package:ensembldb':
filter, selectThe following object is masked from 'package:AnnotationDbi':
selectThe following object is masked from 'package:Biobase':
combineThe following objects are masked from 'package:GenomicRanges':
intersect, setdiff, unionThe following object is masked from 'package:GenomeInfoDb':
intersectThe following objects are masked from 'package:IRanges':
collapse, desc, intersect, setdiff, slice, unionThe following objects are masked from 'package:S4Vectors':
first, intersect, rename, setdiff, setequal, unionThe following objects are masked from 'package:BiocGenerics':
combine, intersect, setdiff, unionThe following object is masked from 'package:biomaRt':
selectThe following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, unionhealthy_epireads <- readRDS("/researchers/caitlin.page/cf_nano/r_output/fixed_healthy_epireads.rds")
luad_epireads <- readRDS("/researchers/caitlin.page/cf_nano/r_output/fixed_luad_epireads.rds")- cpg and genom positions on luad samples were wrong
- was a bug in the for loop
- has now been fixed
cg_sites <- cbind(data.frame(Biostrings::matchPattern("CG", BSgenome.Hsapiens.UCSC.hg38[["chr22"]])), seqnames = "chr22") %>%
mutate(pos = paste0(seqnames, "-", start)) %>%
relocate(pos, seqnames)Warning in .local(x, row.names, optional, ...): 'optional' argument was ignoredcg_sites$index <- 1:nrow(cg_sites)read lengths between cancer and normal
columns <- c("read_id", "seqnames", "start", "end", "read_length", "type", "sample")names(healthy_epireads)[1] "hu5.10" "hu5.12" "ispro.bc2" "ispro.bc3" "ispro.bc4" "ispro.bc5"
[7] "hu5.11" names(luad_epireads)[1] "ispro.bc1" "ispro.bc10" "ispro.bc11" "ispro.bc8" "ispro.bc9"
[6] "ispro.s1" read_lengths <- rbind(healthy_epireads$hu5.10 %>% mutate(type = "healthy", sample = "hu5.10") %>% .[,columns],
healthy_epireads$hu5.11 %>% mutate(type = "healthy", sample = "hu5.11") %>% .[,columns],
healthy_epireads$hu5.12 %>% mutate(type = "healthy", sample = "hu5.12") %>% .[,columns],
healthy_epireads$ispro.bc2 %>% mutate(type = "healthy", sample = "ispro.bc2") %>% .[,columns],
healthy_epireads$ispro.bc3 %>% mutate(type = "healthy", sample = "ispro.bc3") %>% .[,columns],
healthy_epireads$ispro.bc4 %>% mutate(type = "healthy", sample = "ispro.bc4") %>% .[,columns],
healthy_epireads$ispro.bc5 %>% mutate(type = "healthy", sample = "ispro.bc5") %>% .[,columns],
luad_epireads$ispro.bc1 %>% mutate(type = "luad", sample = "ispro.bc1") %>% .[,columns],
luad_epireads$ispro.bc10 %>% mutate(type = "luad", sample = "ispro.bc10") %>% .[,columns],
luad_epireads$ispro.bc11 %>% mutate(type = "luad", sample = "ispro.bc11") %>% .[,columns],
luad_epireads$ispro.bc8 %>% mutate(type = "luad", sample = "ispro.bc8") %>% .[,columns],
luad_epireads$ispro.bc9 %>% mutate(type = "luad", sample = "ispro.bc9") %>% .[,columns],
luad_epireads$ispro.s1 %>% mutate(type = "luad", sample = "ispro.s1") %>% .[,columns])
head(read_lengths) read_id seqnames start end read_length
1 7f095597-d2c5-47a2-93aa-d84b7f9abc94 chr22 10521100 10521455 355
2 e9abf215-8192-48f4-a50f-ef92241e5de3 chr22 10522473 10522658 185
3 4203dd9d-9bff-4d14-b046-7716315e9699 chr22 10524645 10524786 141
4 baf45276-b37b-43fd-8a33-6727de6255fa chr22 10527072 10527210 138
5 a75ce5a8-42f2-4fe2-b62b-35aba109507d chr22 10529269 10529453 184
6 b5e822b7-e149-4ae2-8937-fe4d570a274e chr22 10534032 10534204 172
type sample
1 healthy hu5.10
2 healthy hu5.10
3 healthy hu5.10
4 healthy hu5.10
5 healthy hu5.10
6 healthy hu5.10read_lengths %>% group_by(type) %>% summarise(n())# A tibble: 2 × 2
type `n()`
<chr> <int>
1 healthy 262627
2 luad 496064summary(read_lengths$read_length) Min. 1st Qu. Median Mean 3rd Qu. Max.
45.0 159.0 172.0 228.9 209.0 42104.0 read_lengths %>%
ggplot(aes(x = read_length, colour = type)) +
geom_density(linewidth = 1.5) +
scale_x_continuous(limits = c(0,500)) +
theme_classic()Warning: Removed 33143 rows containing non-finite outside the scale range
(`stat_density()`).
read_lengths %>%
ggplot(aes(x = log2(read_length), colour = type)) +
geom_density()
read_lengths %>%
ggplot(aes(x = log2(read_length), colour = sample)) +
geom_density() +
facet_wrap(~type)
read_lengths %>%
ggplot(aes(x = read_length, colour = type)) +
geom_density() +
scale_x_continuous(limits = c(0,500))Warning: Removed 33143 rows containing non-finite outside the scale range
(`stat_density()`).
* dineika and max expected to see bigger shift to smaller read lengths
in luad but nice to see it in the dinuc * also if you cut with cartesian
coordinates
plot alpha across a region
bc11
healthy_epireads$hu5.10[110:125,] read_id chrom alignment_start alignment_end
110 13570327-b332-4220-b354-a146e96727a9 chr22 10718212 10718857
111 9876388b-0232-4727-bd7e-a552842a33da chr22 10718486 10718651
112 8577ed6e-958e-40a4-9b9c-9ae7013ed69e chr22 10718510 10718832
113 82187fab-97d0-4b38-a392-748c7a21ed4e chr22 10718851 10719039
114 d41898fa-3f6c-40bd-b3ec-00f26ba8509b chr22 10719262 10719407
115 4e0cd750-9d21-46c8-bcce-36b8094e17c6 chr22 10719280 10719402
116 2ac67824-31c7-426f-94a7-a414da518353 chr22 10719305 10719721
117 7608d674-97be-4a7e-804e-fb65de516f13 chr22 10719944 10720514
118 4c9e1d7a-81b2-4f1e-a579-100cff5933cc chr22 10720335 10720501
119 a406a0cd-eb07-46a9-b1ef-75017200a52d chr22 10720590 10720872
120 c7820a2a-d59e-443e-bd0c-109f5f1e8022 chr22 10720773 10720922
121 64ec736a-28f2-49da-a185-5f734791638d chr22 10720790 10721100
122 38c4b052-13d1-400d-8633-485662183d58 chr22 10721104 10721416
123 5f43b35f-e1ae-4f83-a622-7ab1f168a0c8 chr22 10721329 10721549
124 b05df0f3-c19f-4f27-b0dc-b341db141d66 chr22 10721628 10721837
125 6a2a3405-4468-48b2-abc9-95903a14f067 chr22 10721733 10721909
read_length num_cg meth_pattern seqnames start end width
110 645 15 TCCTTTCCCCCCCTC chr22 10718212 10718857 646
111 165 4 TTCC chr22 10718486 10718651 166
112 322 11 TCTTTTCCTCC chr22 10718510 10718832 323
113 188 4 TCCT chr22 10718851 10719039 189
114 145 5 CCCCC chr22 10719262 10719407 146
115 122 4 CCTT chr22 10719280 10719402 123
116 416 12 CCTTCTCCTTCT chr22 10719305 10719721 417
117 570 13 TCCCCCTTTTTTT chr22 10719944 10720514 571
118 166 2 TT chr22 10720335 10720501 167
119 282 4 TTTT chr22 10720590 10720872 283
120 149 1 C chr22 10720773 10720922 150
121 310 3 TTT chr22 10720790 10721100 311
122 312 3 TTT chr22 10721104 10721416 313
123 220 3 CCT chr22 10721329 10721549 221
124 209 7 TCCTTTC chr22 10721628 10721837 210
125 176 5 TTCTC chr22 10721733 10721909 177
index_cg
110 2045
111 2049
112 2049
113 2060
114 2069
115 2070
116 2071
117 2089
118 2100
119 2104
120 2108
121 2108
122 2111
123 2114
124 2118
125 2121
cpg_positions
110 2045,2046,2047,2048,2049,2050,2051,2052,2053,2054,2055,2056,2057,2058,2059
111 2049,2050,2051,2052
112 2049,2050,2051,2052,2053,2054,2055,2056,2057,2058,2059
113 2060,2061,2062,2063
114 2069,2070,2071,2072,2073
115 2070,2071,2072,2073
116 2071,2072,2073,2074,2075,2076,2077,2078,2079,2080,2081,2082
117 2089,2090,2091,2092,2093,2094,2095,2096,2097,2098,2099,2100,2101
118 2100,2101
119 2104,2105,2106,2107
120 2108
121 2108,2109,2110
122 2111,2112,2113
123 2114,2115,2116
124 2118,2119,2120,2121,2122,2123,2124
125 2121,2122,2123,2124,2125
genom_positions
110 10718297,10718327,10718377,10718407,10718577,10718583,10718602,10718637,10718673,10718687,10718707,10718727,10718737,10718752,10718798
111 10718577,10718583,10718602,10718637
112 10718577,10718583,10718602,10718637,10718673,10718687,10718707,10718727,10718737,10718752,10718798
113 10718867,10718889,10718917,10718977
114 10719277,10719302,10719327,10719357,10719382
115 10719302,10719327,10719357,10719382
116 10719327,10719357,10719382,10719407,10719433,10719442,10719483,10719517,10719577,10719607,10719632,10719693
117 10719963,10720037,10720072,10720142,10720162,10720182,10720207,10720217,10720277,10720283,10720307,10720472,10720497
118 10720472,10720497
119 10720602,10720652,10720692,10720697
120 10720892
121 10720892,10720922,10720996
122 10721151,10721211,10721216
123 10721416,10721446,10721521
124 10721656,10721681,10721707,10721742,10721746,10721771,10721821
125 10721742,10721746,10721771,10721821,10721846
num_meth total alpha
110 10 15 0.6666667
111 2 4 0.5000000
112 5 11 0.4545455
113 2 4 0.5000000
114 5 5 1.0000000
115 2 4 0.5000000
116 6 12 0.5000000
117 5 13 0.3846154
118 0 2 0.0000000
119 0 4 0.0000000
120 1 1 1.0000000
121 0 3 0.0000000
122 0 3 0.0000000
123 2 3 0.6666667
124 3 7 0.4285714
125 2 5 0.4000000healthy_epireads$hu5.10 %>% filter(start >= 10718212, end <= 10721909) %>%
ggplot(aes(x = start, y = alpha, colour = read_id)) +
geom_segment(aes(x = start, xend = end, y = alpha, yend = alpha)) +
theme(legend.position = "none")
luad_epireads$ispro.bc11 %>% filter(start >= 10718212, end <= 10721909) %>%
ggplot(aes(x = start, y = alpha, colour = read_id)) +
geom_segment(aes(x = start, xend = end, y = alpha, yend = alpha)) +
theme(legend.position = "none")
* nice that there’s a difference between the cancer and healthy * bc11
has highest tumour purity
maybe beta?
- I don’t really know what I should be doing…
- I think dineika needs team meeting let’s set some goals
- I don’t want to try and run deconvolution stuff
- other people can haggle with scripts and then I just take their instructions and run with it
- which does kinda leave the more exploratory stuff again
- or I could bin the genome - though if jovana was going to do that?
let’s just combine all the samples but have more columns
healthy_epireads[3]